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1.0 

INTRODUCTION 


An  FRC  formation  experiment  could  provide  a  test  bed  for  ionization  kinetics  models  since 
neutral  species  can  propagate  easily  across  field  lines,  but  ionized  species  can  not.  Thus,  the 
depth  of  penetration  of  the  ionized  species,  as  identified  by  spectral  methods,  is  a  useful 
diagnostic  for  determining  the  validity  of  such  models. 

This  report  describes  work  by  NumerEx  to  simulate  the  formation  of  an  FRC  using  mach2.  The 
goal  was  only  to  demonstrate  the  capability  of  doing  the  simulations,  not  to  simulate  any  actual 
experimental  configuration.  Some  very  general  details  of  a  possible  experimental  configuration 
were  inferred  from  the  proposal:  “Magnetized  Target  Fusion  -  A  Proof-of-principle  Research 
Proposal”  by  K.F.  Shoenberg  and  R.E.  Siemon  [LA-UR-98-2413].  This  proposal  describes  an 
experimental  program  to  investigate  the  feasibility  of  Magnetized  Target  Fusion  (MTF)  by  liner 
implosion  of  an  FRC  target.  The  first  phase  of  this  project,  the  formation  of  the  FRC,  will  be 
carried  out  at  Los  Alamos  National  Laboratory.  mach2  is  the  primary  tool  used  for  MHD 
simulation  at  AFRL;  however,  the  traditional  kinds  of  problems  for  which  MACH2  is  well  adapted 
to  solving  do  not  include  inductively-driven  formation  schemes  of  the  type  that  will  be  used  for 
the  FRC.  NumerEx  has  completed  work  on  adapting  the  code  to  solve  this  new  class  of 
problems.  Section  II  will  briefly  describe  the  modifications  that  were  made  to  MACH2  and 
explain  how  its  capabilities  have  been  extended.  Section  III  will  show  some  preliminary  results 
from  simulations  of  FRC  formation  based  on  a  sample  configuration  created  to  roughly  match 
the  examples  shown  in  the  above  mentioned  proposal. 


2.0 

CODE  MODIFICATIONS 

It  has  been  difficult  to  simulate  configurations  where  the  plasma  current  is  driven  inductively  by 
time-varying  externally  applied  vacuum  fields  with  mach2.  The  problem  manifests  itself  in  the 
vacuum  region  between  the  plasma  and  the  applied  field.  If  the  vacuum  region  is  modeled  as 
high  resistivity  plasma,  the  magnetic  field  solver  may  diverge  there,  but  if  the  resistivity  is  not 
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set  high  enough  then  the  electric  field  cannot  penetrate  into  the  plasma  where  it  is  needed  to 
drive  the  current. 

The  problem  was  solved  by  rewriting  the  diffusion  solver  to  use  the  flux  function,  t/r  =  r\ ,  in 
place  of  the  in-plane  components  of  the  magnetic  field.  These  field  components  are  then  defined 
in  terms  of  the  flux  function  by  taking  the  curl: 

B  =  VxA0  =S7y/xd/r  (0.1) 

with  the  unit  vector  in  the  out-of-plane  direction.  Hence,  the  field  diffusion  equation 

—  =  -Vxt7VxB  (0.2) 

dt 

is  replaced  by  the  flux  diffusion  equation 

^  =  Jir2V-r-2V^.  (0.3) 

dt 

This  has  two  important  advantages.  First,  since  the  resistivity  only  appears  outside  the 
differential  operator,  it  can  be  divided  out,  leading  to 

G^L  =  r^-r-^W  (0.4) 

dt 

Since  only  the  conductivity  need  be  specified,  it  can  be  set  equal  to  zero  in  the  vacuum  region, 
making  the  resistivity  infinite  there.  This  is  the  equation  that  the  new  code  solves  with  a  fully 
implicit  Euler  step.  The  conjugate  gradient  method  is  used  to  solve  the  implicit  equations.  The 
multi-grid  solvers  used  elsewhere  in  the  code  might  be  more  efficient  but  they  are  too  complex  to 
have  been  rewritten  in  the  time  allowed  for  this  project. 

The  second  important  advantage  of  using  the  flux  function  is  that  the  magnetic  field  is 
automatically  divergence  free  because  it  is  defined  as  the  curl  of  a  vector  potential.  This 
property  can  be  reproduced  in  the  numerical  equations  if  the  finite  volume  difference  operators 
satisfy  div(curl)=0.  The  original  operators  in  MACH2  did  not  satisfy  this  exactly,  but  it  was  found 
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that  they  could  be  made  to  do  so,  by  changing  the  geometry  coefficients  from  which  they  are 
computed.  The  version  of  the  code  used  for  the  simulations  described  here  maintains  div(B)=0 
exactly.  An  explanation  of  how  this  is  done  will  be  written  up  and  submitted  for  publication  at  a 
later  time. 

The  Lagrangian  hydro  step  in  mach2  does  not  have  to  be  altered  because  that  step  must  conserve 
flux  on  physical  grounds;  hence,  the  flux  function  does  not  change  across  it.  This  observation 
makes  it  possible  to  retain  the  implicit  hydro  step  and  its  stabilizing  iteration  on  the  magnetic 
field,  yet  discard  the  resulting  magnetic  field  in  favor  of  the  flux  function  at  the  beginning  of  the 
iteration.  The  magnetic  field  advection  is  then  replaced  by  advection  of  the  flux  function  itself. 

MACH2  now  has  the  capability  of  simulating  problems  with  any  kind  of  time-varying  vacuum 
field  and  can  do  pure  vacuum  field  diffusion  calculations  of  the  type  that  were  formerly  done 
using  the  D  code1. 


3.0 

SIMULATION  RESULTS 

Magnetic  Fields  The  formation  scheme  simulation  uses  a  20  cm  long,  6  cm  radius  main  field 
coil  and  two  5  cm  long,  6  cm  radius  cusp  coils  at  each  end  of  the  main  coil,  but  separated  from  it 
by  5  cm  long  insulating  gaps.  The  configuration  is  symmetric  about  the  mid-plane  and  only  the 
left  half  is  shown  in  all  the  plots.  The  simulation  included  the  complete  device  and  both  halves 
remained  symmetric.  The  flux  in  each  coil  is  specified  as  a  function  of  time  to  be  a  linear  ramp 
up  to  3  /ts  and  constant  thereafter.  This  corresponds  to  a  square  wave  voltage  pulse  lasting  3  /is. 
Of  course,  any  waveform  could  have  been  specified,  or  the  external  circuit  could  have  been 
specified  and  the  voltage  computed  self-consistently.  The  flux  scheme  does  not  permit  direct 
specification  of  the  current  in  the  coils.  Rather,  it  must  be  determined  from  the  simulation  by 
integration.  The  flux  values  are  set  so  that  the  initial  reversed  bias  field  is  -4.7  kG  and  it  swings 
to  +33  kG  in  3  /is  (with  no  FRC  present— it  is  higher  with  the  FRC  because  the  same  flux  is 
compressed  into  a  smaller  area).  The  cusp  field  starts  at  5.6  kG  and  swings  to  40  kG  in  the  same 
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3  ix s.  The  cusp  field  is  always  slightly  higher  than  the  main  field  to  confine  the  FRC.  Figure  1 
shows  the  coil  geometry  and  the  initial  bias  field.  The  dimensions  are  in  meters.  The  right  edge 
of  the  plot  is  the  mid-plane  and  the  left  edge  is  one  end.  Symmetry  boundary  conditions  are  used 
for  the  field  at  the  ends. 
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Figure  1.  Flux  Plot  Showing  Initial  Bias  Field  at  t  =  0. 

Physical  Model  The  physical  model  consists  of  magnetic  diffusion,  magnetohydrodynamics, 
anisotropic  thermal  diffusion,  and  radiation  cooling  by  emission.  Mass  is  allowed  to  flow  out 
through  the  ends  and  the  single-temperature  energy  equation  includes  Ohmic  heating,  radiation 
losses  and  anisotropic  thermal  conduction  with  thermal  boundaries  treated  as  heat  sinks  held  at  1 
eV.  The  initial  mass  density  is  2xl0'5kg/m3  of  deuterium  for  most  of  the  runs  described  here. 

This  corresponds  to  an  initial  ion  number  density  of  6.7xl015cm  3 .  One  run  uses  half  this  value 
and  another  uses  a  quarter.  The  initial  temperature  is  always  1  eV. 

Anomalous  Restivity  Spitzer  resistivity  is  never  adequate  for  this  type  of  simulation  because  it 
invariably  leads  to  hot,  highly  conducting  vacuum  regions  which  prevent  electric  fields  from 
diffusing  into  the  plasma.  Some  type  of  anomalous  resistivity  model  is  needed  to  prevent  low 
density  regions  from  becoming  highly  conductive.  Unfortunately,  the  state  of  the  art  in 
anomalous  resistivity  theory  has  not  developed  to  the  point  where  any  particular  model  is  known 
to  be  applicable  to  a  given  situation.  The  only  guidance  available  is  to  look  for  past  experiments 
where  some  empirical  information  may  have  been  acquired.  Past  FRC  experiments  at  LANL 
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have  shown  that  the  Chodura  model  with  empirically  determined  parameters  has  led  to  some 
agreement  with  simulation2  so  this  model,  with  their  parameters,  was  used  for  this  simulation. 
Three  runs  were  done  using  this  model,  with  three  different  initial  densities.  One  run  was  also 
done  using  the  lower  hybrid  model  with  MACH2’s  default  parameters  and  one  run  was  done  with 
Spitzer,  but  with  the  Ohmic  heating  turned  off  in  regions  with  density  less  then  5xl014cm’3. 

Discussion  and  Conclusions  Figures  2  through  1 1  are  flux,  and  number  density  contour  plots 
showing  the  time  evolution  of  an  FRC  formation  simulation  using  the  Chodura  model  with  initial 
density  ni=  6.7xl015. 


The  contour  line  that  goes  to  r  =  0  is  <J)  =  0,  and  the  region  it  encloses  to  the  lower  right  has 
negative  flux  values.  Note  that  the  flux  contour  interval  is  much  smaller  for  negative  values  than 
for  positive.  Other  simulations  with  the  different  anomalous  resistivity  models  and  the  different 
initial  densities  look  very  similar  with  only  slight  differences  in  the  size  of  the  FRC. 
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Figure  2.  Flux  and  Density  Plots  at  t  =  0.5  fxs  for  the  Chodura  Model  with  m  =  6.7xl015. 
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Figure  3.  Flux  and  Density  Plots  at  t  =  1.0  ns  for  the  Chodura  Model  with  nj  =  6.7xl015. 
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Figure  4. 


Flux  and  Density  Plots  at  t  =  1.5  ns  for  the  Chodura  Model  with  n;  =  6.7xl015. 
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Figure  5.  Flux  and  DensityPlots  at  t  =  2.0  ns  for  the  Chodura  Model  with  n;  =  6.7xl015. 
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Figure  6.  Flux  and  Density  Plots  at  t  =  2.5  \xs  for  the  Chodura  Model  with  ni  =  6.7xl015. 
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Figure  7.  Flux  and  Density  Plots  at  t  =  3.0  j is  for  the  Chodura  Model  with  ri[  =  6.7xl015. 
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Figure  8.  Flux  and  Density  Plots  at  t  =  3.5  (is  for  the  Chodura  Model  with  n(  =  6.7xl015. 
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Figure  9 . 


Flux  and  Density  Plots  at  t  =  4.0  ps  for  the  Chodura  Model  with  nj  =  6.7xl015. 
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Figure  10. 


Flux  and  Density  Plots  at  t  =  4.5  jxs  for  the  Chodura  Model  with  m  =  6.7xl015. 
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Figure  11. 


Flux  and  Density  Plots  at  t  =  5.0  jis  for  the  Chodura  Model  with  =  6.7xl015. 


There  are  significant  differences,  however,  in  the  peak  densities  and  in  the  temperatures.  Figures 
12,  13  and  14  show  results  from  the  three  different  anomalous  resistivity  models  at  3  fxs. 
Compare  the  Ohmic-cutoff  model  in  Figure  12  to  the  Chodura  in  Figure  14  and  notice  the 
smaller  size  and  60%  higher  density.  Figure  13  shows  the  lower  hybrid  model  at  3  [is;  compare 
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this  to  Figure  14  and  notice  the  intermediate  size  and  density.  The  dependence  of  the  density  on 
the  choice  of  anomalous  resistivity  model  is  most  clearly  seen  in  Figure  15  which  shows  a  time 
history  of  the  peak  density  in  the  three  different  models.  Note  the  differences  in  the  equilibrium 
density  levels  after  3  jtts  and  the  collapse  in  density  after  8  /as  for  the  Chodura  model  only.  These 
variations  limit  the  confidence  that  can  be  placed  on  any  density  prediction. 

No  Anomalous  Resistivity.  No  Ohmic  Heating  below  5x1 O20. 

3.0  microseconds 


0.06 


z 


Figure  12.  Flux  and  Density  Plots  at  t  =  3.0  jas  for  the  Ohmic-cutoff  model  with  nj  = 
6.7xl015. 
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Figure  13.  Flux  and  Density  Plots  at  t  =  3.0  jas  for  the  Lower  Hybrid  Model  with 
nI=  6.7xl015. 
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Figure  14.  Flux  and  Density  Plots  at  t  =  3.0  |js  for  the  Chodura  Model  with  n;  =  6.7xl0ls. 
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Figure  15.  Peak  density  Time  Histories  for  Three  Anomalous  Resistivity  Models,  All  with 
nr=  6.7xl015. 

Figures  16  and  17  show  results  at  3  /x s  from  runs  using  the  Chodura  model  with  different  initial 
density.  They  look  very  similar  except  for  a  slightly  smaller  size  for  the  lower  densities.  Figure 
18  shows  the  time  history  of  the  peak  density  in  these  three  cases.  The  ratios  of  the  equilibrium 
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densities  follow  the  ratios  of  the  initial  values  but  the  density  collapse  occurs  much  sooner  at  the 
lower  densities.  The  physics  of  this  density  collapse  has  not  been  investigated. 


z 


Figure  16.  Flux  and  Density  Plots  at  t  =  3.0  fis  for  the  Chodura  Model  with  n\  =  3.3xl015. 
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Figure  17.  Flux  and  Density  Plots  at  t  =  3.0  j is  for  the  Chodura  Model  with  ii[  =  1.7xl015. 
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Maximum  (on  Density 
Three  Different  Initial  Densities 
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Figure  18.  Peak  Density  Time  Histories  for  the  Chodura  Model,  nj  =  6.7xl015, 3.3xl015, 
1.7X1015. 

Figures  19,  20  and  21  show  temperature  plots  for  each  of  the  three  different  anomalous 
resistivity  models,  all  with  the  same  initial  density.  They  are  all  about  the  same  in  the  FRC 
interior,  ranging  from  around  200  eV  for  the  Chodura  and  lower  hybrid  models,  down  to  150  eV 
for  the  Ohmic-cutoff  model. 
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No  Anomalous  Resistivity.  No  Ohmic  Heating  below  5  x  1 020. 
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Figure  19.  Temperature  Plot  at  t  =  3.0  |xs  for  the  Ohmic-Cutoff  Model  with  ni  =  6.7xl015. 
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Figure  20.  Temperature  Plot  at  t  =  3.0  ps  for  the  Lower  Hybrid  Model  with  m  =  6.7xl015. 
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Figure  21.  Temperature  Plot  at  t  =  3.0  ps  for  the  Chodura  Model  with  nt  =  6.7xl015. 

Figures  21,  22  and  23  show  temperature  plots  for  the  Chodura  model  at  three  different  initial 
densities.  The  temperature  of  the  FRC  increases  by  about  the  same  ratio  as  the  density  drops: 
200  eV  to  400  eV  to  800  eV.  This  makes  the  pressure  about  the  same  for  each  case,  which  it 
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must  be  because  in  equilibrium  it  is  balanced  by  the  magnetic  field  pressure,  which  is  the  same. 
The  hottest  FRC  decays  the  fastest  -  because  of  the  density  dependence  in  the  anomalous 
resistivity  —  if  these  results  can  be  believed.  The  question  of  whether  this  anomalous  resistivity 
will  be  applicable  to  any  experiment  cannot  be  answered  until  the  experiment  is  carried  out. 
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Figure  22.  Temperature  Plot  at  t  =  3.0  ps  for  the  Chodura  Model  with  n(  =  3.3xl0is. 
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Figure  23.  Temperature  Plot  at  t  =  3.0  ps  for  the  Chodura  Model  with  ni  =  1.7xl015. 
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4.0 

CONCLUSION 

Simulations  of  FRC  formation  experiments  can  now  be  performed  using  the  2*/2-dimensional 
magnetohydrodynamic  code  MACH2.  This  report  describes  the  modifications  to  the  code 
employed  and  simulations  performed  to  test  the  model  that  explore  three  different  anomalous 
resistivity  models  and  three  different  initial  densities  for  one  of  those  models.  The  maximum 
density  and  temperature  depend  weakly  on  the  resistivity  model,  and  strongly  on  the  initial 
density. 
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